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Abstract 

Background: Haematotoxicity of conventional chemotherapies often results in delays 
of treatment or reduction of chemotherapy dose. To ameliorate these side-effects, 
patients are routinely treated with blood transfusions or haematopoietic growth 
factors such as erythropoietin (EPO) or granulocyte colony-stimulating factor (G-CSF). 
For the latter ones, pharmaceutical derivatives are available, which differ in absorption 
kinetics, pharmacokinetic and -dynamic properties. Due to the complex interaction of 
cytotoxic effects of chemotherapy and the stimulating effects of different growth factor 
derivatives, optimal treatment is a non-trivial task. In the past, we developed 
mathematical models of thrombopoiesis, granulopoiesis and erythropoiesis under 
chemotherapy and growth-factor applications which can be used to perform clinically 
relevant predictions regarding the feasibility of chemotherapy schedules and 
cytopenia prophylaxis with haematopoietic growth factors. However, interactions of 
lineages and growth-factors were ignored so far. 

Results: To close this gap, we constructed a hybrid model of human granulopoiesis 
and erythropoiesis under conventional chemotherapy, G-CSF and EPO applications. 
This was achieved by combining our single lineage models of human erythropoiesis 
and granulopoiesis with a common stem cell model. G-CSF effects on erythropoiesis 
were also implemented. Pharmacodynamic models are based on ordinary differential 
equations describing proliferation and maturation of haematopoietic cells. The system 
is regulated by feedback loops partly mediated by endogenous and exogenous EPO 
and G-CSF. Chemotherapy is modelled by depletion of cells. Unknown model 
parameters were determined by fitting the model predictions to time series data of 
blood counts and cytokine profiles. Data were extracted from literature or received 
from cooperating clinical study groups. Our model explains dynamics of mature blood 
cells and cytokines after growth-factor applications in healthy volunteers. Moreover, we 
modelled 1 5 different chemotherapeutic drugs by estimating their bone marrow 
toxicity. Taking into account different growth-factor schedules, this adds up to 33 
different chemotherapy regimens explained by the model. 

Conclusions: We conclude that we established a comprehensive biomathematical 
model to explain the dynamics of granulopoiesis and erythropoiesis under combined 
chemotherapy, G-CSF, and EPO applications. We demonstrate how it can be used to 
make predictions regarding haematotoxicity of yet untested chemotherapy and 
growth-factor schedules. 
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Background 

Haematotoxicity of conventional multidrug, multi-cycle chemotherapies often results 
in delays of treatment or reduction of chemotherapy dose [1]. To ameliorate these 
side-effects, patients are routinely treated with blood transfusions (erythrocyte, platelet 
concentrates) or haematopoietic growth factors. For the latter ones, highly potential 
pharmaceutical derivatives are available, which however differ in their pharmacokinetic 
and -dynamic (PK/PD) properties. Due to the complex interaction of cytotoxic effects of 
chemotherapy and the stimulating effects of different growth factors, optimal treatment 
is a non-trivial task. 

In the past, we developed mathematical models of thrombopoiesis (formation of 
platelets responsible for coagulation), granulopoiesis (formation of neutrophil granulo- 
cytes responsible for unspecific immune defense) and erythropoiesis (formation of red 
blood cells responsible for oxygen supply) under chemotherapy and growth-factor appli- 
cations which can be used to perform clinically relevant predictions regarding feasibility 
of chemotherapy schedules and cytopenia prophylaxis with haematopoietic growth fac- 
tors [2-4]. A few predictions were tested and validated in subsequent clinical trials, 
resulting in improved granulocyte colony-stimulating factor (G-CSF) schedules [5,6]. 
However, interactions of lineages and growth- factors were ignored so far. 

In this paper, we aim to combine our models of erythropoiesis and granulopoiesis by 
constructing a model of both cell lines under chemotherapy, G-CSF and erythropoi- 
etin (EPO) applications. Rather than constructing a model from scratch we rely on the 
above mentioned established models to perform a second step towards a more com- 
prehensive model. This implies that we kept the major model assumptions, equations 
and parameters of the isolated erythropoiesis and granulopoiesis model as far as 
possible. 

Our goal is to create a hybrid model of human granulopoiesis and erythropoiesis 
usable for simulations of various chemotherapy schedules with EPO and G-CSF 
support taking into account interactions of the lineages via G-CSF effects on ery- 
thropoiesis [7-9]. Combining the models at stem cell level is a challenging issue, 
since the regulation of the stem cell compartment differs between the single lineage 
models. 

In the present paper, we describe adaptations necessary to combine the models in detail. 
We compare the new hybrid model with the single lineage models for a number of sce- 
narios such as cell loss or single injections of growth-factors. We also validate the hybrid 
model on the basis of a large number of data sets obtained from phase I, II and III clin- 
ical trials containing growth-factor or chemotherapy applications or combinations of it. 
Finally, we perform predictions regarding combined EPO and G-CSF applications during 
chemotherapy. 

Methods 

General structure of the model 

We aim at constructing a hybrid model of human erythropoiesis and granulopoiesis 
under chemotherapy, G-CSF and EPO applications by combining an ordinary differential 
equations model of erythropoiesis [10] and granulopoiesis [3] established recently. The 
single lineage models have a similar structure: Both consist of a number of concatenated 
cell compartments representing different proliferating and maturing cell stages in bone 
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marrow and circulation (see Figure 1). Dynamics of cell compartements are described by 
balance equations of the general form: 

-C = A^^.d^- - (1) 
dt T 

where A^^ is the amplification of the influx derived from the overall amplification A (see 
Additional file 1 section A.l for details). T is the transition time, is the efflux rate 
from the preceding cell compartment, and C the content of the cell compartment (see 
[3,4,10-14]). 

The quantities A and T are usually regulated by growth factors, namely EPO for 
the erythropoiesis model and G-CSF for the granulopoiesis model. Concentrations of 
these growth-factors as well as external applications are also explicitely modelled. The 
quantities A and T are regulated between a minimum and a maximum according to the 
following sigmoidal functions {Ye {A, T} in the following) 

y = Zy (CcytoCO; Y^'^ r--, r--, by) (2) 



I — Scellflux 
^= cytokine flux 
— ► regulation 
inhibition, cell 
destruction 

Oxygen Diss. Curve 
Atmospheric O2 




Figure 1 Structure of the model of erythropoiesis and granulopoiesis under chemotherapy, G-CSF 
and EPO application. Model compartments are presented in boxes (S = stem cells, BE = burst forming 
units - erythroid, CE = colony forming units - erythroid, PEB = proliferating erythrocytic blasts, MEB = 
maturing erythrocytic blasts, RET = reticulocytes, ERY = erythrocytes, EPO = erythropoietin, CG = 
granulopoietic progenitor cells (colony forming units of granulocytes and macrophages), PGB = proliferating 
granulopoietic precursor cells (myeloblasts, promyelocytes, myelocytes), MGB = maturing granulopoietic 
precursor cells (metamyelocytes, banded and segmented granulocytes)), GRA = granulocytes, LY = 
lymphocytes, CX = chemotherapy. Several regulatory feedback loops are displayed. The most important two 
are mediated by EPO and by G-CSF which are produced endogenously and could also be applied externally. 
Chemotherapy is modelled by a transient depletion of cells. 
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(3) 



where Y^'^^ is the corresponding steady state value, Ccyto ^ {Qp^o' ^GCSf} 
centration of either EPO or G-CSF and by is the sensitivity of Y under stimulation. Y^^^ 
and Y^^^ are minimum and maximum values of Y respectively. The function 3 is called 
regulatory function in the following (see also [12], p. 69). 

Amplification and maturation time may depend on the concentration of the growth 
factors G-CSF and EPO, denoted as C^^^pCO for the G-CSF concentration in the central 
compartment, or C^^q(0 for the EPO internalised by red blood cells respectively. Details 
of growth factor mediated regulations as well as assumptions regarding pharmacokinetics 
and -dynamics of growth-factors are explained in sections 'Erythropoiesis model' and 
'Granulopoiesis model'. 

Both single lineage models contain similar but not identical models of stem cell dynam- 
ics which have a different structure than Equation 1. Since the stem cell compartment 
is crucial for combining the models, we explain it in more detail in section 'Stem cell 
compartment S'. Finally, a model of chemotherapy effects was introduced to the lineage 
models which is explained in section 'Chemotherapy model'. A complete set of all model 
equations can be found in the Additional file 1. 

Erythropoiesis model 

The cell kinetic model of erythropoiesis is adopted from earlier modelling works of 
our group with respect to erythropoiesis in mice and humans [10-13,15,16]. This part 
of the model describes the development of mature erythrocytes from haematopoietic 
stem cells and its regulation by endogenous EPO. It consists of the cell compartments 
S (stem cells), BE (burst forming units - erythroid), CE (colony forming units - ery- 
throid), PEB (proliferating erythrocytic blasts), MEB (maturing erythrocytic blasts), RET 
(reticulocytes), FRY (erythrocytes), and EPO concentrations of different sites (see also 
Figure 1). 

EPO is assumed to increase the proliferation and to shorten the maturation time of red 
blood cells in the bone marrow [12,14,16]. Endogenous production of EPO (EPOprod) [16] 
is assumed to depend on the oxygen partial pressure in the kidneys and the number of cir- 
culating red blood cells [11,13]. A detailed description of the model of EPO applications 
can be found in [10]. In brief, it consists of a pharmacokinetic model of EPO (adapted 
from [17]) and a subcutaneous injection model of EPO (adapted from [18]). The PK 
model consists of three compartments: EPO in central serum (Qpq), the peripheral com- 
partment (^C^pQ^, and EPO bound to receptors (C^pq) [17], with first order transitions 
between the peripheral and the central compartment describing reversible protein bind- 
ing of EPO. EPO can bind to free EPO receptors forming a drug-receptor complex which 
dissociates again or EPO is internalised. EPO is removed from circulation by the latter 
mechanism or via unspecific elimination modelled by a first order loss term. Internalised 
EPO serves as argument of all regulation functions regarding EPO. 

Different pharmaceutical EPO derivatives are modelled: EPO Alfa, EPO Beta, EPO 
Delta and Darbepoetin Alfa. The first three are recombinant human EPO derivatives of 
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the first generation differing only in their glycosylation patterns. Darbepoetin Alfa is a 
next generation EPO pharmaceutical with increased half-life in serum. 

Different application sites of EPO result in different absorption kinetics [18]. While 
intravenous injections can simply be modelled by pulse functions, subcutaneous injec- 
tions are characterised by delayed and incomplete absorption resulting in a reduced 
bioavailability [19-21]. To account for this fact, we adapted a model of subcutaneous EPO 
injection developed for sheep [18]. The model contains two absorption routes: directly 
from the subcutaneous tissue into the bloodstream, or indirectly via the lymphatic sys- 
tem. In both processes, we included time delays by concatenated sub-compartments with 
first order transitions (see Additional file 1, section A. 3 for details). A loss of EPO is 
assumed at the injection site and in the lymphatic system. The structure of this injec- 
tion model is described in detail in [10,18]. All equations can be found in the Additional 
file 1. 

Due to differences in the lymphatic flow at different anatomical regions [18,19,22], we 
split available injection modi into different groups. While the parameters describing sub- 
cutaneuos injection can differ between these groups, the PK parameters are the same for 
endogenous EPO and all EPO derivatives except for Darbepoetin. Available data allowed 
to specify parameter settings for injections of EPO Alfa into thigh, EPO Alfa into shoul- 
der, EPO Alfa into forearm, EPO Alfa into upper arm or abdomen, EPO Beta into forearm, 
EPO Beta into arm or abdomen, EPO Beta into thigh, EPO Delta, and finally, Darbepoetin 
Alfa (see [10]). 

A model of chemotherapy action on erythropoiesis is added as explained in 
section 'Chemotherapy model'. 

Granulopoiesis model 

The cell kinetic model of granulopoiesis also constructed by our group [3,4,12,23] 
describes the development of mature granulocytes from haematopoietic stem cells and 
its regulation by G-CSE It consists of the cell compartments S (stem cells), the granu- 
lopoietic progenitors CG (colony forming units of granulocytes and macrophages), the 
proliferating granulopoietic precursors PGB (myeloblasts, promyelocytes, myelocytes), 
the maturing granulopoietic precursors MGB (metamyelocytes, banded and segmented 
granulocytes), the granulocytes GRA and G-CSF concentrations at different sites (see also 
Figure 1). 

We consider both, endogenously produced G-CSF and injections of the derivatives Fil- 
grastim and Pegfilgrastim. Endogenous G-CSF production is regulated by the demand of 
mature granulocytes. In analogy to former versions of our model [3,4], this is modelled 
phenomenologically as a function of the content of the final bone marrow compartment 
and circulating granulocytes. 

Pharmacokinetics of endogenous G-CSF and injections are modelled by three com- 
partments: a subcutaneous compartment C^ ^^^ in which G-CSF pharmaceuticals are 
injected, a central compartment C^^^p in which the drugs are pharmacologically active 
and a peripheral compartment C^^^^p representing reversible protein binding [24]. This 
model was originally developed for mice and rats and later adopted for the human 
situation [3,25,26]. 

Exogenous G-CSF applications are modelled by an injection function analogous to EPO 
injections. Delayed influx of injected G-CSF from the subcutaneous compartment into 
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the central compartment caused, e.g., by lymphatic absorption [18] is modelled by divi- 
sion of the subcutaneous compartment into two subcompartments (for details regarding 
modelling of delays see Additional file 1 section A.3). Dose-dependent bioavailability of 
G-CSF [26] is modelled by a Michaelis-Menten loss term within the first sub compartment 
of the subcutaneous tissue [3]. Transitions between central and peripheral compartments 
were modelled by two-way first order kinetics [27] . 

G-CSF is eliminated from the central compartment via two routes: Unspecific renal 
elimination of G-CSF is modelled by a first order kinetic [27,28]. Specific degradation 
mediated by the number of circulating granulocytes is modelled by a Michaelis-Menten 
kinetic proportional to this number [3,25,29-32] . 

We use the same pharmacokinetic and -dynamic parameter settings for Filgras- 
tim and endogenous G-CSF because of their high similarity [33-35]. The observed 
differences between Filgrastim and Pegfilgrastim [4,28,36-40] are modelled using the 
same model structure but different parameters settings for absorption, distribution, 
degradation and regulatory mechanisms (regulatory functions). Receptor competi- 
tions between Pegfilgrastim and endogenous G-CSF (or Filgrastim) is modelled by 
adding the regulatory functions of Pegfilgrastim and Filgrastim using a weighting 
factor. 

Pharmacokinetic and -dynamic modelling of G-CSF is described in detail in [3,4]. In 
our former model, we assumed a delayed effect of G-CSF on bone marrow cells. Since this 
effect was rather small, we decided to drop it in the combined model in order to reduce 
complexity and computational burden. 

A model of chemotherapy is attached to our granulopoiesis model as explained in 
section 'Chemotherapy model'. 

Stem cell compartment S 

Cells differentiating into granulopoietic or erythropoietic lineages originate from a com- 
mon pool of cells called haematopoietic stem cells ([41]). We used the same equations 
of stem cell dynamics in our models of granulopoiesis and erythropoiesis, but regula- 
tions were slightly different. Basic concepts can be traced back to the book of Loeffler & 
Wichmann [12]. Since the stem cell compartment is crucial for combining the models we 
present its structure in more detail now. 

The stem cell compartment has self-renewal capability. To achieve a steady state, on 
average 50% of the proliferating stem cells remain in this compartment, and 50% differen- 
tiate into red or white blood cell lineages. Hence, the compartment equation differs from 
Equation 1: 



d ac. 
at rs 



/^OUt 



2{l-p)Cs 



as 



(4) 
(5) 



where rs is the average duration of a cell cycle, p the self-renewal probability, as is the 
proliferative fraction, Cs is the size of the stem cell compartment and C^^* its efflux com- 
mitting to the haematopoietic lineages. According to [12], the self-renewal probability p 



Schirm et al. Theoretical Biology and Medical Modelling 201 4, 1 1 :24 
http://www.tbiomed.eom/content/1 1 /I /24 



Page 7 of 28 



is regulated by the demand of the hematopoietic bone marrow system. In steady state, we 
assume p^^^ = ^. Thus, for the initial conditions (steady state) it holds that 

Cs(0) = Cr = 1 (6) 

q^t^O) = Cg'^^-^''' = 2 (1 -/?^°') (7) 

Similar to the stem cell model of [12], self-renewal of stem cells p is regulated by a 
competition of the stem cell content (C^\t)), and the granulopoietic (C^^\t)) or the 
erythropoietic (^C^\t)) bone marrow cells. 

p=p (cf(t). cf(t). c^\t).P8. ^E, ^G, m)) , 

The parameters i^e>^g> and are hypothetical weighting factors originally defined 
as the strength of the influence of the bone marrow cells C^\t) = C^\t) = 

Cbe(0+Qe(^)+Qeb(0+C'meb(0 ^rel/^\ _ Cqg (0+Qgb(0+C'mgb 

/<-'nor_i_/<-'nor_i_/<-'nor _i_/<-'nor > dllU. J — i<^nor_i_i<^nor lienor • 

^BE ^^CE ^^PEB^^BE ^ '^CG ~'~^PGB~'~^MGB 

According to [12], it is assumed that 

= ^nor _ ^min ^ ^max _ ^nor 



m) = 



rel/^\0.6 



cf(ty 



for Cf(t) 
for Cf(t) 



< 1 
> 1 



p=P8 tanh (-i^sit) (cf(t) (cfit) -l)-i^G {C^\t) - l)) + 0.5, 

(8) 

where p""""' = ^, = -8, t^e = -2. 

In our single lineage models, this regulation was simplified in different ways. Since ery- 
thropoiesis was not included in the model of granulopoiesis, only granulopoietic cells C^^ 
were assumed to influence the self-renewal probability in S with the factor t^q = —10, i.e. 
we assumed C^^ = C^^ and thus 

p=P8 tanh(-?^s(0 (Cf(t) - l) - i^G {C'^\t) - l)) + 0.5. 

In contrast, a constant granulocyte value Cq^ = 1 was assumed in our erythropoiesis 
model. The erythropoietic bone marrow cells C^^ were multiplied by i^e = —2, thus 

p=P8 tanh (-m) (cfit) - l) - (cfit) - l)) + 0.5. 

The proliferative fraction can be interpreted as the percentage of cells which are 
currently in cell cycle. The proliferative fractions ax of the compartments S, BE or CG are 
also regulated by the haematopoietic bone marrow system O^'^it) and €^^{1) or 0^^(t): 

^ ^ I /^rel//.\ ^rel/^\ ^rel/^\ ^min ^nor _int ^max \ 
(^X = '^^X I v^^' yClx ,Clx yClx , <^E> <^G> I > 
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Again, the parameters cos.coq, and co^ are weighting factors. They represent the strengths 
of the influence of stem cells, erythropoietic and granulopoietic cells on the proliferative 
fraction of BE, CG and S respectively. With 

x = coEln Cf (t) + COG In C^^^ (t) (9) 



InCfit), for q'^<l 
Cf(t) - 1, for cf > 1 



y= \ — In ^ - In 



the proliferative fraction is given by 



ax 



e-y+ey ^ ^X ^ ^X ^ ^ ^X QO) 
for at^ 



■^x — X — X — X 



It is a monotone function ranging between a^^"^ and a^^^. Low cell numbers in the 
bone marrow compartments cause a higher demand of proliferating cells, and therefore, 
a larger proliferative fraction ax- The value of y defines the actual point on the regulatory 
curve. The variable ;v is a measure of the total bone marrow content. It is calculated as a 
weighted sum of the logarithms of the relative counts of stem cells, erythropoietic cells 
and granulopoietic cells. If any cell counts tend to zero, x tends to minus infinity, and 
with it, a becomes maximal. Parameter values a}"^^ corresponds to x = — In 2 and a^^^ 
corresponds to = 0 (see Figure S2 in the Additional file 1). 

In analogy to the regulation of the self-renewal probability, a factor coq = 0.4 is 
assumed for the influence of the granulopoietic bone marrow compartments C^^ in the 
granulopoiesis model, thus 



\nCf(t), for Cf <1 
Cf(t) - 1, for cf > 1 



In contrast, in the erythropoiesis model we assumed a constant granulocyte value 
Cq^ = 1 and the logarithmized erythropoietic bone marrow cells Cf are multiplied by 
0)^ — 0.3, i.e. 

1 ^rel. ^ \ lnCf(tl for Cf < 1 ] 

E V ^ ^ [ Cf(t) - 1, for cf > 1 J 
Chemotherapy model 

We used the same model of chemotherapy in both the erythropoiesis and the granu- 
lopoiesis model. It is based on the following assumptions which are extensively discussed 
in [2-4]: Chemotherapy results in a delayed, reversible and transient depletion of the 
compartments S, CG, PGB, MGB, BE, CE, FEB, MEB, and RET, where the toxic effect 
is specific for the cell stages, for different drugs or drug combinations, and for different 
doses of the same drug. This is quantified by corresponding sets of toxicity parameters 
(see tables A.8 and A.9 in the Additional file 1). A higher toxicity in the first chemotherapy 
cycle is modelled by a factor f^^^^ > 1. It is assumed that different cytotoxic drugs damage 
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:Cx = .C^-^-^f^-Cx (11) 



independently of each other. This allows us to add toxicity functions of chemotherapeutic 
drugs applied in combinations. 

Infusion of chemotherapeutic drugs is again modelled using pulse functions. The effect 
of chemotherapy is introduced to the balance equations of the bone marrow cell com- 
partments by a first-order loss term ^x- Hence, the modified schematic compartment 
Equation 1 has the form 

-Cx=A 

for X e {CG, PGB, MGB, BE, CE, PEB, MEB, RET}, ^^x depends on the kind, dose and 
timing schedule of chemotherapy as well as on the affected bone marrow comparted X 
(see Additional file 1 for further details). Analogously, the modified stem cell Equation 4 
has the form 

^Cs = (2;7-l)Cs-->I's-Cs. (12) 
at rs 

In clinical practice often only leukocytes, i.e. the sum of unspecific (granulocytes) 

and specific (lymphocytes) immune cells, are available. To avoid additional modelling of 

lymphopoiesis, the reduced lymphocyte count under chemotherapy is described by an 

exponential function of the form (see [3]): 

Cgra (t) 



CwBC (t) ^ Cix exp (-^LY (0) + CGRA- 



i^nor 
^GRA 



where cly = 3000 cells per /x/ and cgra = 4000 cells per /x/ are the steady state con- 
centrations of lymphocytes and granulocytes respectively, ^ly is the toxicity function for 
lymphocytes which is analogously defined as the toxicity functions of bone marrow cell 
stages. 

Age is a major risk factor of haematotoxicity under chemotherapy [1]. In our for- 
mer granulopoiesis model [4], we assumed different toxicity parameters for the age 
groups >60 years and <60 years. This assumption will be carried over to our hybrid 
model. 



Combination of granulopoiesis and erythropoiesis models 

It is our major intention to keep the single lineage models as unchanged as possible. But 
a few assumptions are necessary to combine them. These assumptions refer to the reg- 
ulation of stem cells, the interaction of granulopoiesis and erythropoiesis, and common 
chemotherapy effects. 

In the following, we describe the assumptions and necessary modifications of the single 
lineage models in detail. A schematic representation comprising all parts of the model is 
shown in Figure 1. 

Common stem cell compartment 

In the single lineage models, the concurrent influence of erythropoietic and myeloid 
cells on the regulation of the stem cell compartment, i.e. the regulation of self-renewal 
probability and proliferative fraction, could not be considered. This resulted in a num- 
ber of simplifications explained in section 'Stem cell compartment S'. When combining 
the single lineage models, these simplifications are no longer required. Hence, we rely 
on the extended regulatory functions regarding proliferative fraction and self-renewal 
(Equations 8, 9 and 10). The same regulatory principles were used for the proliferative 
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fractions of the compartments BE and CG, i.e. these quantities now also depend on both, 
relative stem cell count (C^^), relative size of bone marrow erythropoiesis (C^^) and 
granulopoiesis (C^^). 

Furthermore, we implemented a splitting of the output of the stem cell compartment 
which feeds both lineages. We assume that 15% of the released cells (C^^^, see Equation 7), 
differentiate into red blood cells = 0.15), and 80% into the white blood cell line 
(ofG = 0.8, see also [12], pp. 61, 62). However, these parameters are dummies having no 
impact on the behaviour of the model in its current form since only relative changes of 
compartment sizes are considered. They will become relevant later when implementing a 
model of lineage commitment [42]. Hence, the influx into the subsequent erythropoietic 
and granulopoietic progenitor compartments reads as follows: 

= o.i5q^^ 

CS, = 0.8Cs-^ 

Note that 5% of stem cells are assumed to commit to thrombopoietic lineage. 
Compartment CE and G-CSF effect on erythropoiesis 

In our single lineage model of erythropoiesis, the amplification and transition time in 
compartment CE is only regulated by the growth factor EPO. Now we additionally assume 
an inhibiting influence of endogenous G-CSF, Filgrastim and Pegfilgrastim on this com- 
partment. Evidence for interactions of Erythropoietin and G-CSF and mutual influence 
of the different lineages is described in [7,8,43]. 

This interaction is modelled by additional regulatory functions of G-CSF multiplied to 
those of amplification and transition time of the compartment CE (see Equation 3 and 
Figure 2). More precisely, we have 

Ace = Zac^ {C^^o) ' ^a_gcsf(0 (13) 

TcE = Ztc^ (C^^p^o) • ^t_gcsf(0. (14) 

with 

i=A_GCSF(0 = Wp(t) • ZA_Peg(0 + (1 - Wp(t)) • ZA_Fil(0 (15) 

^T_GCSf(0 = Wp(t) • ZT_Peg(0 + (1 - Wp(t)) . Zj_m(t) (16) 

where Zac^> Zt^^ are the regulatory functions of amplification and transition time in 
compartment CE^ Zt Fii(0> Za Fil(0 ^rid Za Peg(0> Zt Peg(0> ^re regulatory functions of 
endogenous G-CSF respectively of Filgrastim and Pegfilgrastim, 0 < wp(t) < 1 is the 
weighting factor to model the superimposing effect of concurrent endogenous G-CSF / 
Filgrastim and Pegfilgrastim action (see [3]), and C^^q is the internalised EPO (see [10]). 
The fitting process performed later resulted in a rather small influence of Filgrastim and 
Pegfilgrastim on the amplification in CE. Therefore, we decided to set Fa_gcsf(0 = 1. A 
complete set of parameters regarding this regulation is listed in Table 1. 

Common chemotherapy model 

The same chemotherapy model was used in our single lineage models of granulopoiesis 
and erythropoiesis. However, combination of the models required some adjustments of 
toxicity parameters due to the interactions of the lineages not considered in the previous 
models. 
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Figure 2 Regulatory functions in compartment CE. Regulatory functions describing tine effect of Filgrastim, 
Pegfilgrastim, EPO and combinations of it on amplification rate and transition time in the compartment CE. 
While amplification is assumed to be unaffected by G-CSF, the transition time is delayed under G-CSF 
stimulation. A: Delay factor of transition time - dependence on endogenous G-CSF/Filgrastim concentration. 
We mark the values achieved at steady-state (green) and at maximum after a single s.c. application of 480 /xg 
(black). B: Delay factor of transition time - dependence on Pegfilgrastim concentration. Values achieved at 
maximum after a single s.c. application of 6000 /xg (black) and steady state (green) are marked. C: 
Amplification in CE - dependence on internalised EPO. No G-CSF effects are assumed here. D: Transition times 
in CE - dependence on internalised EPO and G-CSF concentration. We present the raw regulatory functions 
without considering G-CSF effects (dotted) and maximal deviations after either Filgrastim application of 480 
fig (solid) or Pegfilgrastim application of 6000 /xg (grey). Note that the regulatory functions presented in A and 
B are superimposed (Equation 1 6) under Pegfilgrastim or combined Pegfilgrastim/Filgrastim applications . 



Numerical methods for simulation 

Simulations were performed with MATLAB 7.5.0.342 (R2007b) using the SIMULINK 
toolbox (The MathWorks Inc., Natick, MA, USA). Numerical solutions of the equation 
system are obtained using the variable step solver from Adams and Bashford (odellS, 
SIMULINK toolbox). 

Data sets 

We combined the data sets used to establish the single lineage models of erythropoiesis 
and granulopoiesis. These data sets comprise time courses of absolute neutrophil counts 



Table 1 Parameters for modelling G-CSF effects on compartment CE 



Parameter 


Meaning 


Filgrastim 




Pegfilgrastim 




T-min T-min 
' CEf ' ' CEp 


Factor of transition time under minimal stimulation 


0.9995 


Fitted 


1 


Set 


-rnor -rnor 
'CEf' 'CEp 


Factor of transition time under normal stimulation 


1 


Set 


1.037 


Fitted 


T-max T-max 
'cEf' ' CEp 


Factor of transition time under maximal stimulation 


98.37 


Fitted 


2.787 


Fitted 


'CEf' 'CEp 


Sensitivity of factor of transition time 


0.930 


Fitted 


0.5660 


Fitted 
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(ANC), white blood cell counts (WBC), hemoglobin (HB), hematocrit (HCT), serum 
concentration of EPO and G-CSF, red blood cell counts (RBC), reticulocyte counts (RET), 
or percentages of reticulocytes after EPO or G-CSF application in healthy volunteers and 
time series data of hemoglobin, ANC or WBC in patients treated with chemotherapy with 
or without G-CSF or EPO support. Data sets without access to raw data were taken from 
literature. The automated tool "ycasd" [44] was used to extract the data as precisely as 
possible. However, in literature often only data of one blood cell lineage per scenario were 
available. 

We also have access to raw patients data of the German High Grade Non-Hodgkin's- 
Lymphoma Study Group [45-48] and the German Hodgkin's Lymphoma Study Group 
[49,50], and the German Breast Group [51]. These data comprise both, haemoglobin and 
leukocyte data for 11 different chemotherapeutic drug combinations. Considering dif- 
ferent growth factor schedules, we analysed raw data of 15 different therapy scenarios. 
Patients HB data were censored after receiving erythrocyte concentrates. For some ther- 
apy schedules we split the data taking into account individual risk groups (young vs. 
elderly patients). Combined with the literature data, we modelled 15 different chemother- 
apeutic drugs or drug combinations. Taking into account differences in growth factor 
schedules results in a total of 33 different therapy schemes. An overview of all data sets 
used for modelling is presented in table A. 10 in the Additional file 1. 

Scenarios were either used to determine unknown parameters of the model or to 
validate model predictions. 

Estimation of parameters 

It is our goal to keep the parameters of our single lineage models constant as far as possible 
to reduce the problem of over-fitting, but a few adaptations were necessary due to the 
changes required for model combination (see section 'Combination of granulopoiesis and 
erythropoiesis models'). While the parameters of the cell kinetic models of granulopoiesis 
and erythropoiesis and the pharmacokinetics of EPO and G-CSF described in [3,4] and 
[10] were kept constant, we modified chemotherapy parameters and those involved in 
regulating amplification and transition time in the modified compartment CE (see tables 
of chemotherapy parameters in the Additional file 1 and Table 1). 

To estimate toxicity parameters, data of at least one scenario are required in which the 
corresponding drug or drug combination was used. Hence, chemotherapy parameters 
were established using the EC-T and E-T-C scenarios with G-CSF and with or without 
EPO [51], CHOP 21 and CHOEP 21 for young and elderly patients [45,46], the BEACOPP 
21 and BEACOPP 21 escalated data [49,50], the high CHOEP 21 data [47], the literature 
data with ESHAP and 12 x 5 /xg/kg Filgrastim or 100 /xg/kg Pegfilgrastim on day 6 [52] 
and the literature data with Platinum and Etoposide with and without Darbepoetin [53]. 
Toxicity parameters for the literature scenarios without HB data (Doxorubicin-Docetaxel 
[54-57] or Carboplatin-Paclitaxel [58]) were taken from the granulopoiesis model. Due 
to lack of data, toxicity parameters for erythropoiesis could not be established for these 
scenarios. 

Parameters regarding Filgrastim effects on the compartment CE were determinated 
using the HB- and WBC-data of the NHL-B trial (CHOP 14, elderly patients. Filgrastim 
480 /xg day 4-13), the data of E-T-C with EPO and G-CSF on day 3-10, and the data of 
the "Ricover" trial (elderly patients, CHOP 14, G-CSF Filgrastim 480 /xg day 6-12), since 
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these scenarios have simple chemotherapy regimens but different G-CSF schedules infor- 
mative for our purposes. Parameters regarding Pegfilgrastim effects on the compartment 
CE were determinated on the basis of the data of CHOP 14 with Pegfilgrastim on day 2 
and day 4 [5]. 

For parameter fitting, (l+3)-evolutionary-strategies with self-adapting mutation step 
size [59,60] were used. This refers to a non-deterministic optimisation algorithm with 
one possibly immortal parent and three children per generation. We applied the fitness 
function described in [3,10] based on the area between simulation and data curve: 



/model is the simulation result at time t for the parameter set k = /q, . . ./c^. to < 
t < ti describe the range of time points for which data are available, and /data (0 is the 
interpolated curve of data medians. Since we model relative sizes of cell compartments, 
model outputs are multiplied by normal values to allow comparisons with clinical data. 

The fitness value is defined as the left hand side of Equation 17. It is specific for the clin- 
ical intervention (e.g. application of a certain chemotherapy), the outcome measure (e.g. 
granulocyte counts) and the clinical data set (e.g. median of patient data at certain time 
points) considered. Fitness values of different outcomes of the same clinical intervention 
are added to calculate an overall fitness value of the scenario. For simultaneous fitting of 
more than one scenario, we add corresponding fitness values. 

We performed a sensitivity analysis of newly introduced cell kinetic parameters of 
our model by changing their values by 2.5% and calculating corresponding deteriora- 
tions of the fitness function. It revealed that minimum and normal values of the G-CSF 
regulations of compartment CE (see Table 1) have considerably higher precision than cor- 
responding estimates of maximum values and b-parameters (details see Additional file 1: 
Figure S4). 

Validation of model 

Most of the scenarios for which we have data were not used for parameter fitting but 
for validation purposes. This requires that necessary parameters for the simulation of a 
particular scenario were either determined by previous fitting steps or can be taken from 
the single lineage models. 

Throughout, scenarios with growth-factor administration in healthy volunteers were 
used to validate our hybrid model (data sets 1-46, 61-84 in table A.IO of the Additional 
file 1), since the majority of required parameters were determined in the framework of 
the single lineage models. 

The situation for scenarios including chemotherapy is more complicated: At least one 
scenario must be fitted to determine the toxicity parameters of an unknown drug combi- 
nation, but often, toxicity parameters determined in the single lineage models performed 
also fine for the hybrid model. This applied for the scenarios 85-96 of table A.IO com- 
prising literature data of Carboplatin and Paclitaxel or Doxorubicin and Docetaxel with 
injections of Filgrastim or Pegfilgrastim or CHOP-like therapies (therapies which are 
similar to the CHOP combination therapy consisting of the drugs Cyclophosphamide, 
Doxorubicin, Vincristin and Prednisone). Additionally, the data of high CHOEP 14 (data 
set 57), CHOEP 14 and CHOP 14 (data sets 53 and 55, subgroup of young patients) were 




(17) 
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also not included into the fitting procedure but used for validation. For the majority of 
chemotherapy scenarios, only WBC or ANC data are available. Hence, validation of the 
granulopoiesis sub-model relies on a considerably larger data base. 

Results 

Model behaviour 

After calibration, we study the qualitative behaviour of the model and compare it later 
with the results of the single lineage models. In our hybrid model, we have additional 
feedbacks due to the G-CSF effect on CE cells, the indirect effect of erythropoietic bone 
marrow cells on granulopoiesis and vice versa via regulation of the proliferative frac- 
tions of S, CE, BE and CG as well as the regulation of stem cell self-renewal probability. 
Feedbacks are negative throughout resulting in more or less damped oscillations of com- 
partment sizes after disturbance of the system. This behaviour is robust against changes of 
model parameters except for those involved in stem cell feedback: If proliferative fraction 
and stem-cell self-renewal are allowed to respond more intensively to changes of bone 
marrow contents, disturbances of the system result in stable oscillations, a phenomenon 
which was formerly used to explain cyclic neutropenia in dogs ([61,62], see Figure S3 in 
the Additional file 1). 

We study model dynamics after single perturbation and compare them with the results 
of the single lineage models (see section 'Comparison of the single lineagemodels and the 
hybrid model'). At first, we simulate a single i.v. injection of 150 lU/kg EPO Alfa at t=0 
(see Figure 3): The MEB compartment of the hybrid model expands quickly, followed by 
a peak in reticulocytes. Since initially, amplification in BE and CE is stimulated by EPO, 



single injection of 150 lU EPO/kg 




0 5 10 15 20 25 30 35 



days 

Figure 3 Model behaviour after single injection of EPO. We present cell counts normalised to steady 
state values after i.v. injection of 1 50 lU/kg EPO Alfa. After damped oscillations of compartment sizes the 
system converges to equilibrium state. 
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these compartments also show increasing cell counts. The compartment S also increases 
initially, caused by the up-regulation of the self-renewal probability due to increased 
bone marrow cell counts. In contrast, the compartment PEB is depleted by the short- 
ened maturation time and thereby increased efflux into the compartment MEB. After 
a minimal decrease in the number of CG and PGB the granulopoietic cells follow the 
oscillation in S. Overall, damped oscillations of the compartment sizes can be observed. 
Next, we simulated a single s.c. injection of Filgrastim at t=0 (see Figure 4 left panels): 
ANC grows fast caused by the increased release of mature granulocytes from bone mar- 
row to circulation. Thereafter, the GRA compartment is reduced rapidly caused by the 
fast degradation of Filgrastim and the short half life of granulocytes. Since the prolifera- 
tion in PGB is increased by G-CSF, this compartment also grows initially. The number of 
stem cells increases (similarly to EPO application) due to up-regulation of self-renewal. 
The number of erythropoietic cells is decreased due to the G-CSF effect on CE. Again, 
all cell counts show damped oscillations, which is typical for short-term disturbances of 
the system. After Pegfilgrastim administration (see Figure 4 right panels), the cell num- 
bers react similarly, but increments and reductions of cell counts remain for a longer time 
frame than after Filgrastim application. In contrast to Filgrastim application, no oscil- 
lations occur after Pegfilgastim application. The reason for this behaviour is the longer 
half-life of Pegfilgrastim compared to Filgrastim. Pegfilgrastim applications result in 
prolonged stimulations of the bone marrow, and with it, an expanded stem cell pool which 
is only gradually reduced after elimination of Pegfilgrastim. 



single injection of 300 |j.g Filgrastim 



single injection of 300 ^g Pegfilgrastim 




Figure 4 Model behaviour after single injection of G-CSF. Model behaviour after perturbation witli 
respect to steady state values. We present cell counts normalised to steady state values after a single s.c. 
injection of 300 /xg Filgrastim (left) or 300 /ig Pegfilgrastim (right). A single Filgrastim injection results in 
damped oscillations of compartment sizes, which was not observed for Pegfilgrastim. 
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Now we study the effects of continuous stimulation with EPO and G-CSF (see Figure 5). 
For this purpose, we keep EPO and G-CSF serum concentrations constant at higher val- 
ues. This is achieved by constant influxes into these compartments. EPO stimulation 
results in rapidly increasing numbers of MEB and RET. After oscillation, a new steady 
state of BE, CE, MEB and RET on a higher level compared to the unperturbed steady 
state is achieved. In contrast, the cell number in PEB reaches a slightly lower equilib- 
rium after an initial strong decline. After an initial peak the stem cells stabilise on a 
slightly lower level. The numbers of CG, PGB, MGB and ANC decline to a lower steady 
state. 

A continous administration of G-CSF results in an increased equilibrium of PGB, MGB, 
ANC and S. Compartments CG, BE, CE, PEB, MEB, RET and FRY stabilise on lower 
levels. 

Finally, we want to study the system behaviour after a single CHOP chemotherapy 
administration at ^ = 0. Figure 6 illustrates the reaction of the system to this damaging 
effect. After an initial growth of the numbers of CE and PEB, damped oscillations of both 
lineages appear over a longer time period. The numbers of BE, MEB and RET show a 
strong decline after chemotherapy administration. Later, they show damped oscillations 
as well while they slowly approach their normal values. Erythrocytes are slightly dimin- 
ished after chemotherapy and later their number grows slowly without oscillations. The 
compartments S, CG, PGB, MGB and ANC were strongly depleted after chemotherapy 
application. Later their numbers oscillate and they slowly converge to their normal values. 



continuous stimulation 1 0 ID EPO/kg per day continuous GCSF stimulation (30 [ig/d) 




04 I \ [ \ \ 1 0 I ^ ^ ^ ' ' ^ 1 1 

0 100 200 300 400 500 0 50 100 150 200 250 300 350 400 

days days 

Figure 5 Model behaviour during continuous stimulation with EPO or G-CSF. Left: EPO, right: Filgrastim. 
New steady states are readied after a certain time. 
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Figure 6 Model behaviour after a single chemotherapy administration. We present cell counts 
normalised to steady state values after a single administration of CHOP chemotherapy. Except for 
erythrocytes, all lineages show damped oscillations over a longer time period. 



Comparison of the single lineage models and the hybrid model 

In this section, we compare simulation results after perturbations of the hybrid model 
with corresponding results of the single lineage models. As in the previous sections, 
we simulate single injections of EPO, Filgrastim, Pegfilgrastim and CHOP chemother- 
apy (see Figure 7) respectively. Newly implemented interactions between erythropoiesis 
and granulopoiesis in the hybrid model resulted in some differences (for further details, 
see section "Model comparison" in A.9 of the Additional file 1). Generally, in the hybrid 
model, stem cell oscillations are less damped than in the single lineage models, especially 
after chemotherapy (CHOP). Interestingly, despite this observation, the simulation results 
for HB and WBC are comparable between the hybrid and the single cell lineage models, 
i.e. the lack of damping in S is compensated at later cell stages. Nevertheless, we predict 
continued oscillations of mature blood cells after chemotherapy, an issue which could be 
verified in patients. 

Chemotherapy scenarios 

Chemotherapy toxicity parameters were re-estimated by fitting the predictions of the 
model to our data sets (see tables of chemotherapy parameters in the Additional file 1 and 
section 'Estimation of parameters'). 

In Figures 8, 9 and 10 we present a selection of chemotherapy scenarios for which 
WBC and HB data are available, i.e. which could not completely be described by our for- 
mer single lineage models. Scenarios comprise CHOP chemotherapy for the treatment of 
aggressive non-Hodgkins lymphoma disease with different dosing and timing schedules 
of Filgrastim (Figure 8) or Pegfilgrastim (Figure 9) and breast cancer therapies supported 
by Filgrastim and EPO (Figure 10). Note that these are multi-cycle therapies in which 
a drug or drug combination is applied multiple times, usually with 14 days or 21 days 
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Figure 7 Comparison of dynamics after CHOP chemotherapy. Comparison of the behaviour of the 
combined and the single lineage models after perturbation with CHOP chemotherapy. We present stem cell 
dynamics, CE, HB and WBC normalised to steady state values. 



intervals (see figure legends and table A, 10 in the Additional file 1 for details). We 
obtained a good agreement of model and data for all scenarios. Simulation results of the 
single cell lineage models and the combined model are comparable throughout. In the 
Additional file 1 we present the results of further chemotherapy scenarios (Figures S27, 
S33-S40, compare with table A.IO). 

Validation 

Scenarios not included into the fitting procedure are suitable for validation of the hybrid 
model. For example, in [64] time series data of reticulocytes, WBC, HB and HK after mul- 
tiple injections of 5000 lU EPO beta into healthy volunteers are presented. This scenario 
was not included into the fitting procedure of both single lineage models as well as the 
hybrid model. For the latter one, we observed a good fit of the data of both cell lineages 
over a longer time period of 15 weeks (see Figure 11). 

Results of the other validation scenarios can be found in the Additional file 1. We 
observed a reasonable fit for all of these scenarios demonstrating the validity of our 
model. 

Model predictions 

BEACOPP escalated is an intense chemotherapy for the treatment of Hodgkins lym- 
phoma resulting in high degrees of toxicity affecting all haematopoietic lineages [49,65]. 
G-CSF but not EPO support is mandatory. In consequence, anaemia due to cumulative 
toxicity to erythropoiesis is common at later cycles of the therapy. Here we simulate the 



Schirm et al. Theoretical Biology and Medical Modelling 201 4, 1 1 :24 
http://www.tbionned.conn/content/l 1 /I /24 



Page 19 of 28 



CHOP 21. age > 60 



CHOP 21, age > 60 




20 40 60 80 100 120 
CHOP 14, age > 60, Filgrastim day 4-13 
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Figure 8 Model simulations of different CHOP chemotherapies with or without Filgrastim support. 

First row: Six cycles of CHOP witli cycle duration of 21 days without Filgrastim, second row: Six cycles of CHOP 
with cycle duration of 14 days with Filgrastim at day 4-13 at each therapy cycle, third row: Eight cycles of 
CHOP with cycle duration of 1 4 days with Filgrastim at day 6-1 2 at each therapy cycle. We present time series 
of HB and WBC simulated with our combined model (black solid line) and compare it with corresponding 
simulation results of the single lineage models (grey solid line). We also compare model results with our 
clinical data (median: circles, first and third quartile: grey dashed line) taken from the trials published in [45,63]. 
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Figure 9 Model simulations of six cycles of CHOP chemotherapy with cycle duration of 1 4 days and 
different Pegfilgrastim support. We present time series of HB and WBC simulated with our combined 
model (black solid line). We compare model results with our clinical data (median: circles, first and third 
quartile: grey dashed line) taken from clinical trials [5]. 
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Figure 10 Hybrid model simulations after chemotherapy. Filgrastim or EPO. We present results of 
different breast cancer therapies. First row: four cycles of the drug combination epirubicine + 
cyclophosphamide followed by four cycles of paclitaxel with cycle duration of 21 days without G-CSF, second 
row: three single drug cycles of epirubicine, paclitaxel, cyclophosphamide applied consecutively with cycle 
duration of 14 days with Filgrastim on cycle days 3-10, third row: The same therapy as described in second 
row but with additional EPO Alfa. We show HB and WBC values of simulation (black line), medians (circle), 
25th and 75th percentiles (grey dashed line) of patients data. While raw data are available for the first two 
scenarios, data of the third scenario are taken from [51]. 



effect of concomitant G-CSF and EPO applications during the course of the therapy. We 
simulate the application of eight cycles BEACOPP-21 escalated chemotherapy with Fil- 
grastim on days 8-15 with or without weekly application of 300 lU/kg Darbepoetin (see 
Figure 12, left). While WBC dynamics are roughly the same, the decrease of HB is delayed 
by EPO possibly resulting in better quality of life of patients. 

We also analysed the effect of weekly injections of 300 lU/kg Darbepoetin during eight 
cycles of CHOP- 14 chemotherapy in elderly patients. Here, Filgrastim was applied on 
days 6-12 during each cycle of therapy. While dynamics of WBC are similar, HB stabilises 
at higher values during treatment with Darbepoetin. 

Discussion 

Conventional cytotoxic chemotherapy plays a major role in cancer therapy. Development 
of intensified regimen improved the outcome of several diseases [45,46,49,66,67] but is 
limited by toxic side effects. A major dose-limiting side effect is general haematotoxicity 
which is routinely treated with growth factors EPO and G-CSF. Different pharmaceu- 
tical derivatives of these factors are available, which differ greatly in pharmacokinetic 
and -dynamic properties. Furthermore, outcome of growth factor treatment depends on 
many factors such as chemotherapy drugs used, drug doses, growth-factor derivatives and 
individual risk factors [68,69]. Due to this variety of variable therapy parameters, iden- 
tification of optimal growth-factor schedules cannot be performed solely on the basis 
of clinical trials. We showed in the past that mathematical models of haematopoiesis 
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under chemotherapy can facilitate the development of optimised and individualised 
growth-factor schedules [6]. 

Efforts to model haematopoiesis under chemotherapy and growth-factor applications 
are considerable [2,4,11-13,16,23,31,32,70-84]. All these attempts are based on mech- 
anistic or semi-mechanistic parameterical models. In view of the complexity of the 
haematopoietic system, it might be worthwhile to consider model-free forecasts of 
haematopoietic outcomes (see discussion in [85,86]). However, to our knowledge, such an 
approach has never been proposed. 

There are only a few attempts to develop multilineage models or hybrids of exist- 
ing model: In [87] a compartmental model of granulopoiesis and thrombopoiesis 
including G-CSF and high dose chemotherapy application is proposed. Feedback mech- 
anisms from thrombopoietic, granulopoietic and stem cells on stem cells and an 
interaction between the G-CSF compartment and the white blood cells are included. 
Erythropoesis was not considered. A stochastic multilineage model only considering 
haematopoietic stem cells and progenitor cells based on the idea of lineage com- 
mitment can be found in [42]. In [88] an ordinary differential equations model of 
granulopoiesis under chemotherapy and Filgrastim injections was combined with a 
model of stem cell plasticity. G-CSF and chemotherapy effects are modelled but ery- 
thropoiesis was not included. In the recent version of our erythropoeisis model, we 
combined a cell kinetic model with EPO injection models, with an EPO pharma- 
cokinetic model and with a chemotherapy model developed for granulopoiesis. But 
granulopoiesis was not considered. In [89,90] several earlier models of haematopoi- 
etic stem stell dynamics and models for the regulation of neutrophils, platelets and 
erythrocytes were combined to a model of the regulation of the hematopoietic 
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Figure 1 1 Example validation scenario. Data of Lundby et al. [64] after multiple injections of EPO Beta. 
Time series of WBC, HB, Reticulocytes and HK are presented. Solid black curve represents simulation results. 
Circles and grey line represent means and standard deviations respectively [64]. 
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Figure 12 Prediction scenarios. We present predictions regarding BEACOPP-escalated and CHOP 14 
(elderly) supported by Filgrastim and weekly Darbepoetin starting on day 0. Simulation results of WBC and 
HB are shown as black curve. Grey curves show predictions of dynamics without additional Darbepoetin 
applications. Median of data is represented by circles, and grey hatched line describe first and third quartile of 
the data [49,63]. 



system describing periodic chronic myelogenous leukemia and cyclical neutrope- 
nia taking G-CSF effects into account. Chemotherapy was not considered in these 
models. 

Here, we propose a comprehensive combined model of human erythropoiesis and gran- 
ulopoiesis under chemotherapy, G-CSF and EPO applications for the first time. For this 
purpose, we combined a former cell kinetic model of erythropoiesis under chemotherapy 
and EPO applications developed by our group [10-13,16] with a cell kinetic model of gran- 
ulopoiesis under chemotherapy and G-CSF applications [3,4,12,23]. These models were 
chosen for their clinically relevant applications and parameter settings established on the 
basis of large clinical datasets. Models are based on an ordinary differential equations 
system describing major cell stages in bone marrow and blood as well as dynamics of 
key cytokines. Rather than constructing a combined model from scratch we rely on these 
established modeling works and adopted their assumptions, equations and parameter set- 
tings as far as possible. Therefore, this approach can be considered as a second level of 
modeling by the combination of established models in order to build more comprehensive 
ones. 

Accordingly, only a few additional assumptions and only a few new parameters were 
introduced in order to combine the model. The most challenging issue was to com- 
bine the models at stem cell level since regulation of stem cells dynamics depend on 
mature cell stages which was modelled differently in our single lineage models. This 
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was solved by substituting the stem cell regulation by a model originally developed 
by Loeffler and Wichmann [12]. The model includes elements of lineage competition 
with respect to stem cell differentiation. Evidences of interactions between cell lineages 
regarding stem cell committment are summarised in [7]. Since model behaviour depends 
heavily on stem cell dynamics, it was by far not clear whether this approach would be 
successful Indeed, this adaptation later required re-estimations of toxicity parameters 
of chemotherapy since stem cell toxicities are highly sensitive with respect to model 
outcome. 

A second assumption was made with respect to suppression of erythropoiesis by G- 
CSF: Previous experiments have shown, that G-CSF administration leads to a reduction of 
CFU-E in the femur of mice [8,91-93]. Bungart et al. measured a decline of bone marrow 
CFU-E in mice to about % of the control value [91]. Bensinger et al. detected a greater 
decrease of hematocrit values in human blood donors receiving rhG-CSF than in donors 
without G-CSF [9]. However, the groups with or without G-CSF were not comparable. In 
view of these observations, we assumed that G-CSF decreases proliferation and prolonges 
transition time in CE [43]. However, our model simulations show that assuming a delay of 
the transition time in dependance on G-CSF dose and derivative is sufficient to describe 
the phenomena observed in human studies. Generally, suppression of erythropoiesis by 
G-CSF becomes apparent only under high serum levels of G-CSF over a longer time scale 
such as after Pegfilgrastim injections. 

In summary, the major achievement of our present modelling efforts is that we suc- 
cessfully combined our single lineage models of erythropoiesis and granulopoiesis under 
chemotherapy by introducing a common stem cell model and an interaction of the lin- 
eages mediated by G-CSF. The model is mechanistic in the sense that it directly describes 
amplification and maturation processes of cell stages, relevant cell fluxes and PK/PD of 
EPO and G-CSF derivatives. However, we also have to acknowledge that our model is 
semi-mechanistic with respect to stem cell feedback and chemotherapy action which 
both are only phenomenologically described. This constitutes possible starting points for 
future model improvements. 

While we rely on the parameter estimates of the single lineage models in general, 
newly introduced parameters based on assumptions necessary to combine the model 
were determined by fitting the predictions of the model to clinical data. Additionally, as 
mentioned earlier, it was also necessary to re-estimate toxicity parameters of chemother- 
apeutic drugs due to the changes required for stem cell modelling. It was our objective 
to construct a model predicting medians of patients, i.e. patients heterogeneity was not 
considered so far. We considered virtually all available time series data of mature blood 
parameters and cytokine levels after application of G-CSF, EPO, chemotherapy and com- 
binations of it in order to parametrise the model or to validate its results. Data were 
retrieved either from the literature or from our own clinical study groups for which we 
have access to raw data of patients. This results in more than 100 scenarios comprising 15 
different chemotherapeutic regimens. 

In order to compare the results of our hybrid model with those of our single lineage 
models, we performed a number of qualitative simulations, namely single applications 
of G-CSF, EPO or chemotherapy. We observed some differences with respect to the 
dynamics of stem cells and earlier bone marrow compartments. In general, correspond- 
ing oscillations in the hybrid model are stronger and less damped compared to the single 
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lineage models. However, this effect is only predicted for earlier cell stages for which no 
data are available. Due to additional regulations by G-CSF, the dynamics at later cell stages 
are almost identical to those of the former version of the model showing less pronounced 
oscillations at stem cell level. A similar observation was made by combining a two- 
dimensional difference equation model of stem cell regulation with our granulopoiesis 
model [88]. Thus, one can conclude that dynamics of blood cells are less informative 
regarding stem cell dynamics, i.e. appropriateness of different stem cell models can hardly 
be decided on the basis of available clinical data. 

With the hybrid model we simulated more than 100 scenarios including growth fac- 
tor applications in healthy volunteers, chemotherapy without growth factor application 
and chemotherapy with supportive G-CSF, EPO or both. Although there is a consider- 
able inhomogeneity between patient groups considered (different age groups, diseases, 
distributions of risk factors) a good agreement between model simulation and data was 
achieved for almost all scenarios using a single parameter set (except for different toxi- 
city parameters used for young and elderly patients). However, one has to acknowledge 
that the dynamics of white blood cells during chemotherapy is richer, and thus, more 
informative regarding appropriateness of model behaviour than the dynamics of the 
erythropoietic lineage. Reticulocyte data are scarcely available while erythrocytes, haema- 
tocrit or haemoglobin show more or less a constant decline in the time course of the 
therapy. 

The majority of available data sets were used for model validation rather than parame- 
ter estimation. This especially applies for data after growth factor application into healthy 
volunteers for which no additional toxicity parameters are required. Additionally, a num- 
ber of chemotherapy scenarios for which toxicity parameters were available from earlier 
modelling steps, were also used for model validation. The good agreement of model and 
data for almost all validation scenarios shows that the model covers a wide range of 
scenarios. Unfortunately, often data of only one of the lineages are available from the lit- 
erature. Data received from our clinical trials are most detailed resulting in close meshed 
time series at a daily scale supported by dozens of single measurements. However, only 
data of leukocyte and haemoglobin are available to support our modelling. It would 
greatly benefit from more detailed blood parameters measured in parallel such as erythro- 
cytes, hematocrit, hemoglobin, reticulocytes, iron status, and lymphocytes, granulocytes 
rather than leukocytes. 

Finally, we demonstrated how the model could be used to make clinically relevant 
predictions regarding the outcome of different growth-factor schedules after chemother- 
apy. This requires that the toxcity parameters of the therapy are available from previous 
fitting steps. Then, the model can be used to simulate and compare alternative growth- 
factor schedules. We predict for example that both, patients treated with CHOP 14 or 
BEACOPP escalated therapy, would benefit from weekly Darbepoetin treatment without 
reducing the effectiveness of concomittant G-CSF treatment. Clinical trials are required 
to validate these model predictions. 

However, we have to acknowledge that the present model only allows median pre- 
dictions while critical time-courses are clinically more relevant and therapy-limiting. 
Although this aspect is not yet covered, there is a clear perspective towards modelling 
individual data either by fitting parameter sets for patient risk groups or by assuming dis- 
tributions of model parameters. Accordingly, we plan to extend our model in the near 
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future and apply it in order to support improvement and individualisation of growth 
factor therapies. 

Conclusions 

We successfully combined established bio-mathematical models of granulopiesis and 
erythropoiesis under chemotherapy and growth-factor applications by introducing a 
common stem cell compartment and lineage interactions. The model explains data of 
about 100 clinical scenarios including 15 different chemotherapies and six growth-factor 
derivatives. We demonstrated how the model can be used to make clinically relevant 
predictions regarding combined G-CSF and EPO treatment during chemotherapy. 
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